--- title: H Parameter keywords: fastai sidebar: home_sidebar summary: "A reader for the binary H parameter map" description: "A reader for the binary H parameter map" nb_path: "05_hpar.ipynb" ---
{% raw %}
{% endraw %} {% raw %}

memmap_binary[source]

memmap_binary(fname)

{% endraw %} {% raw %}

read_hpar_binary[source]

read_hpar_binary(fname=Path('/luna4/maye/h_parameter_maps/hpar_global128ppd_v1c.bin'))

Reading Paul Hayne's binary global 128 ppd H-parameter map.

Parameters
----------
fname: str, pathlib.Path
    Path to binary 128 ppd map file binary (ending in .bin)
{% endraw %} {% raw %}
{% endraw %} {% raw %}
hpar_path
PosixPath('/luna4/maye/h_parameter_maps/hpar_global128ppd_v1c.bin')
{% endraw %}

Size in bytes:

{% raw %}
size = hpar_path.stat().st_size
size
3303270404
{% endraw %}

Scaling down 80 lats to 70:

{% raw %}
n_lats = 20480 / 80 * 70
n_lats
17920.0
{% endraw %} {% raw %}
npix = 46080 * int(n_lats)
npix
825753600
{% endraw %}

multiply by 32 bits, divide by 8 to get bytes:

{% raw %}
expected_bytes = npix * 32 / 8
expected_bytes
3303014400.0
{% endraw %}

Ratio:

{% raw %}
size / expected_bytes
1.0000775061713325
{% endraw %} {% raw %}
npix = 46081 * (int(n_lats) + 1)
npix
825817601
{% endraw %} {% raw %}
expected_bytes = npix * 32 / 8
size / expected_bytes
1.0
{% endraw %} {% raw %}
lats = memmap_binary(lats_path)
lats[0], lats[-1]
(memmap([69.99609, 69.99609, 69.99609, ..., 69.99609, 69.99609, 69.99609],
        dtype=float32),
 memmap([-70.00391, -70.00391, -70.00391, ..., -70.00391, -70.00391,
         -70.00391], dtype=float32))
{% endraw %} {% raw %}
lons = memmap_binary(lons_path)
lons[:, 0], lons[:, -1]
(memmap([0.003906, 0.003906, 0.003906, ..., 0.003906, 0.003906, 0.003906],
        dtype=float32),
 memmap([360.0039, 360.0039, 360.0039, ..., 360.0039, 360.0039, 360.0039],
        dtype=float32))
{% endraw %} {% raw %}
with np.printoptions(precision=10):
    print(lons[0])
[3.9059999e-03 1.1719000e-02 1.9531000e-02 ... 3.5998828e+02 3.5999609e+02
 3.6000391e+02]
{% endraw %} {% raw %}
with np.printoptions(precision=10):
    print(lons[:, 0])
[0.003906 0.003906 0.003906 ... 0.003906 0.003906 0.003906]
{% endraw %} {% raw %}
np.printopt
{% endraw %} {% raw %}
with np.printoptions(precision=20):
    print(lons[0])
[3.9059999e-03 1.1719000e-02 1.9531000e-02 ... 3.5998828e+02 3.5999609e+02
 3.6000391e+02]
{% endraw %} {% raw %}
np.linspace(-70, 70, int(n_lats))
array([-70.        , -69.99218706, -69.98437413, ...,  69.98437413,
        69.99218706,  70.        ])
{% endraw %} {% raw %}
half_pixel_degree = 1 / 128 / 2
half_pixel_degree
0.00390625
{% endraw %} {% raw %}
with_n_plus_1 = np.linspace(
    -70 - half_pixel_degree, 70 - half_pixel_degree, int(n_lats) + 1
)
{% endraw %} {% raw %}
from_data = lats[
    :,
    0,
][::-1]
{% endraw %} {% raw %}
n = len(from_data)
{% endraw %} {% raw %}
%matplotlib widget
{% endraw %} {% raw %}
from_data[:10]
memmap([-70.00391, -69.99609, -69.98828, -69.98047, -69.97266, -69.96484,
        -69.95703, -69.94922, -69.94141, -69.93359], dtype=float32)
{% endraw %} {% raw %}
with_n_plus_1
array([-70.00390625, -69.99609375, -69.98828125, ...,  69.98046875,
        69.98828125,  69.99609375])
{% endraw %} {% raw %}

class HReader[source]

HReader(lat_limit=None, other_path=None)

Data reader class for H parameter map.

It accesses preproduced data for 128 ppd for DEM, slope, and aspect,
located on the luna4 disk.

It uses virtual dask.Arrays so that virtually no memory is consumed until you
resolve a chain of operations with the `.compute()` call.

Attributes
----------
H: xarray.DataArray
    H parameter map
{% endraw %} {% raw %}
{% endraw %} {% raw %}
H = HReader()
{% endraw %} {% raw %}
H.img
<xarray.DataArray 'H' (lat: 17921, lon: 46081)>
dask.array<where, shape=(17921, 46081), dtype=float32, chunksize=(5792, 5792), chunktype=numpy.ndarray>
Coordinates:
  * lat      (lat) float64 70.0 69.99 69.98 69.98 ... -69.98 -69.98 -69.99 -70.0
  * lon      (lon) float64 0.0 0.007812 0.01562 0.02344 ... 360.0 360.0 360.0
Attributes:
    long_name:  H Parameter
    units:      H units :-P
{% endraw %} {% raw %}
H.get_H_by_coord(20, 50)
0.07237300276756287
{% endraw %} {% raw %}
H.img.sel(lat=slice(21, 20), lon=slice(120,121)).hvplot()
{% endraw %} {% raw %}
H.convert_to_lon180()
{% endraw %} {% raw %}
H.img
<xarray.DataArray 'H' (lat: 17921, lon: 46081)>
dask.array<getitem, shape=(17921, 46081), dtype=float32, chunksize=(5792, 5792), chunktype=numpy.ndarray>
Coordinates:
  * lat      (lat) float64 70.0 69.99 69.98 69.98 ... -69.98 -69.98 -69.99 -70.0
  * lon      (lon) float64 -180.0 -180.0 -180.0 -180.0 ... 180.0 180.0 180.0
Attributes:
    long_name:  H Parameter
    units:      H units :-P
{% endraw %} {% raw %}
H.get_H_by_coord(20, 70)
---------------------------------------------------------------------------
InvalidIndexError                         Traceback (most recent call last)
/tmp/ipykernel_30172/2079116474.py in <module>
----> 1 H.get_H_by_coord(20, 70)

/tmp/ipykernel_30172/1965268766.py in get_H_by_coord(self, lat, lon)
     76     def get_H_by_coord(self, lat, lon):
     77         "Note: Decided to not apply offset!"
---> 78         val = self.img.sel(lat=lat, lon=lon, method="nearest")
     79         return float(val)
     80 

/luna4/maye/miniconda3/envs/py38/lib/python3.8/site-packages/xarray/core/dataarray.py in sel(self, indexers, method, tolerance, drop, **indexers_kwargs)
   1313         Dimensions without coordinates: points
   1314         """
-> 1315         ds = self._to_temp_dataset().sel(
   1316             indexers=indexers,
   1317             drop=drop,

/luna4/maye/miniconda3/envs/py38/lib/python3.8/site-packages/xarray/core/dataset.py in sel(self, indexers, method, tolerance, drop, **indexers_kwargs)
   2472         """
   2473         indexers = either_dict_or_kwargs(indexers, indexers_kwargs, "sel")
-> 2474         pos_indexers, new_indexes = remap_label_indexers(
   2475             self, indexers=indexers, method=method, tolerance=tolerance
   2476         )

/luna4/maye/miniconda3/envs/py38/lib/python3.8/site-packages/xarray/core/coordinates.py in remap_label_indexers(obj, indexers, method, tolerance, **indexers_kwargs)
    419     }
    420 
--> 421     pos_indexers, new_indexes = indexing.remap_label_indexers(
    422         obj, v_indexers, method=method, tolerance=tolerance
    423     )

/luna4/maye/miniconda3/envs/py38/lib/python3.8/site-packages/xarray/core/indexing.py in remap_label_indexers(data_obj, indexers, method, tolerance)
    115     for dim, index in indexes.items():
    116         labels = grouped_indexers[dim]
--> 117         idxr, new_idx = index.query(labels, method=method, tolerance=tolerance)
    118         pos_indexers[dim] = idxr
    119         if new_idx is not None:

/luna4/maye/miniconda3/envs/py38/lib/python3.8/site-packages/xarray/core/indexes.py in query(self, labels, method, tolerance)
    222                     indexer = index.get_loc(label_value)
    223                 else:
--> 224                     indexer = index.get_loc(
    225                         label_value, method=method, tolerance=tolerance
    226                     )

/luna4/maye/miniconda3/envs/py38/lib/python3.8/site-packages/pandas/core/indexes/base.py in get_loc(self, key, method, tolerance)
   3369             tolerance = self._convert_tolerance(tolerance, np.asarray(key))
   3370 
-> 3371         indexer = self.get_indexer([key], method=method, tolerance=tolerance)
   3372         if indexer.ndim > 1 or indexer.size > 1:
   3373             raise TypeError("get_loc requires scalar valued input")

/luna4/maye/miniconda3/envs/py38/lib/python3.8/site-packages/pandas/core/indexes/base.py in get_indexer(self, target, method, limit, tolerance)
   3440 
   3441         if not self._index_as_unique:
-> 3442             raise InvalidIndexError(self._requires_unique_msg)
   3443 
   3444         if not self._should_compare(target) and not is_interval_dtype(self.dtype):

InvalidIndexError: Reindexing only valid with uniquely valued Index objects
{% endraw %} {% raw %}
lons = memmap_binary(lons_path)
{% endraw %} {% raw %}
lons.shape
(17921, 46081)
{% endraw %} {% raw %}
pd.Series(lons[:, 0]).value_counts()
0.003906    17921
dtype: int64
{% endraw %} {% raw %}
pd.Series(lons[:, -1]).value_counts()
360.003906    17921
dtype: int64
{% endraw %} {% raw %}
H = memmap_binary(hpar_path)
{% endraw %} {% raw %}
H.shape
(17921, 46081)
{% endraw %} {% raw %}
H[:, 0][:10]
memmap([0.078986, 0.078389, 0.07891 , 0.078347, 0.07683 , 0.0854  ,
        0.077884, 0.081484, 0.079931, 0.075148], dtype=float32)
{% endraw %} {% raw %}
H[:, -1][:10]
memmap([0.078389, 0.07891 , 0.078347, 0.07683 , 0.0854  , 0.077884,
        0.081484, 0.079931, 0.075148, 0.073657], dtype=float32)
{% endraw %} {% raw %}
H = read_hpar_binary()
{% endraw %} {% raw %}
H.isel(lon=0).hvplot() * H.isel(lon=-1).hvplot()
{% endraw %} {% raw %}
(H.isel(lon=0)[1:] - H.isel(lon=-1)[:-1]).compute().data
array([-0.000521,  0.000563,  0.001517, ...,  0.000944,       nan,
             nan], dtype=float32)
{% endraw %} {% raw %}
left = H.isel(lon=-1)[:-1].compute().data
{% endraw %} {% raw %}
right = H.isel(lon=0)[1:].compute().data
{% endraw %} {% raw %}
pd.Series(left - right).value_counts()
 0.000000    17326
-0.000407        2
-0.000293        2
 0.000192        1
 0.000830        1
             ...  
-0.001465        1
-0.000870        1
-0.001488        1
-0.001567        1
 0.002393        1
Length: 142, dtype: int64
{% endraw %} {% raw %}
17326 / left.shape[0]
0.9668526785714285
{% endraw %} {% raw %}
%matplotlib widget
{% endraw %} {% raw %}
fig, ax = plt.subplots()
ax.plot(H[:, 0], label="left")
ax.plot(H[:, -1], label="right")
[<matplotlib.lines.Line2D at 0x7fed8f00b760>]
{% endraw %} {% raw %}
H[:, 0][:10]
memmap([0.078986, 0.078389, 0.07891 , 0.078347, 0.07683 , 0.0854  ,
        0.077884, 0.081484, 0.079931, 0.075148], dtype=float32)
{% endraw %} {% raw %}
H[:, -1][:10]
memmap([0.078389, 0.07891 , 0.078347, 0.07683 , 0.0854  , 0.077884,
        0.081484, 0.079931, 0.075148, 0.073657], dtype=float32)
{% endraw %} {% raw %}
H.lon.diff("lon") / 2
<xarray.DataArray 'lon' (lon: 46080)>
array([0.00390625, 0.00390625, 0.00390625, ..., 0.00390625, 0.00390625,
       0.00390625])
Coordinates:
  * lon      (lon) float64 0.007812 0.01562 0.02344 ... 360.0 360.0 360.0
{% endraw %} {% raw %}
H = read_hpar_binary()
{% endraw %} {% raw %}
import hvplot.xarray
{% endraw %} {% raw %}
H.isel(lon=0).hvplot() * H.isel(lon=-1).hvplot()
{% endraw %} {% raw %}
rate = 47
{% endraw %} {% raw %}
rate * 8
376
{% endraw %} {% raw %}
arr_x = []
arr_y = []
for i in range(5):
    arr_x.append(np.sort(np.random.uniform(size=20)))
    arr_y.append(np.sort(np.random.uniform(size=20)))
{% endraw %} {% raw %}
plt.figure()
for x, y in zip(arr_x, arr_y):
    plt.plot(x, y)
{% endraw %}